1

Re-use your analysis and run an interaction model between your IV and the other covariate. Plot the interaction as a coefficient plot.
For the coefficient plot, you can stick to easystats.
linear_model_interaction <-
  lm(
    curfew_yes_no ~ age_cat * education_cat,
    data = gp_covid
  )

model_parameters(linear_model_interaction) %>% 
  plot()

2

The plot might be quite vivid and colorful. This one’s rather an open task: Try to make it more dim and more greyish or apply just another theme.
Remember, the plots are base on ggplot2. You could choose from one of ggplot2’s built-in themes here: https://ggplot2.tidyverse.org/reference/ggtheme.html. In another step, you could adjust the scale’s colour, e.g., with scale_colour_grey(). Be creative.
model_parameters(linear_model_interaction) %>% 
  plot() +
  scale_colour_grey(start = 0, end = .5, guide = "none") +
  theme_minimal()
## Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.

3

Plot the interaction effect as a prediction.
For this purpose, you need the plot_model() function from sjPlot.
library(sjPlot)

plot_model(
  linear_model_interaction,
  type = "int"
)